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Abstract 

Correlated systems at both zero and nonzero temperature are treated here from a 
novel angle using a functional method. This functional method is an extension of the 
usual effective potential method. Here, however the effective action is made to depend 
explicitly on the correlation effects that are inherent in the physics involved. This 
will enable us to obtain new expressions for the free energy, the specific heat and the 
ground state energy. The new expansion is shown to give the expected results for the 
homogeneous case at zero temperature. However at non-zero temperature we are able 
to get new sets of diagrams that have a vanishing effect at zero temperature. To lowest 
order these diagrams if summed properly will solve a In T anomaly in the specific heat 
of an electron gas at low temperature. We are also able to show that this method 
provides a very clear way to extend the RPA approximation. We calculate the effect 
of exchange on the ring diagrams at zero temperature and show how to include some 
of the ladder diagrams. Our results agree well with known numerical calculations. We 
conclude by showing that this method is in fact a variant of the time dependent density 
functional method and can in principle be applied to study the effects of correlation in 
the nonhomogeneous case. 
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1 Introduction 



Strongly correlated systems continue to be a popular subject due to the importance of these 
systems in physical processes in general. Electron systems form an important subset of 
these systems. The nonlinearity that appears in this type of problem complicates greatly 
the calculation of any physical properties of the system that strongly depend on the cor- 
relation between the individual particles. In this paper we are mainly concerned with a 
particular aspect of this area, that is, of calculating the ground state equilibrium energy of 
a large Fermi gas in an external potential, at both zero and at nonzero temperature. At 
zero temperature the random phase approximation (RPA) is a well known method which 
addresses this problem. Here we present a way to improve on it, using an effective action 
method, which includes higher order interactions in a consistent way. This new method is 
also applicable at finite temperature. At low temperature, the new method shows clearly 
the importance of the so called "anomalous diagrams" in cancelling the non-physical InT 
term in the specific heat Another important aspect of this new approach is that it is 
applicable to nonhomogeneous Fermi systems. The currently popular method for treating 
such systems is the density functional method (DFT)0, ||. In DFT, the density of the 
system plays a central role and the energy is a functional of it. Both the current method and 
density functional method are extensions of the Thomas-Fermi method. In DFT the density 
is found by solving self-consistently a one-particle Schrodinger equation. Correlation effects 
are taken into account by using the homogeneous result locally. The kinetic energy term is 
taken to be that of a non-interacting system having the same density as the real system under 
consideration. In spite of all these approximations, density functional methods have proved 
to give better values for binding energies than Hartree-Fock in atoms, molecules and solids 
in general. However the method works best only for almost homogeneous systems; surface 
effects are handled poorly by this method. So far the attempt to include a gradient of the 
density in the energy functional makes the result worse. The success of the method is there- 
fore limited. Here we propose a different method, that in principle can treat nonhomogeneity 
in a better fashion. However, the problem becomes harder to manage numerically. 

The method proposed here is a path integral method known as the effective action 
method. In one dimension it reduces to the well-known WKB method to first approxi- 
mation. This work is a first step toward finding a new practical scheme for the evaluation 
of correlation energies and other physical properties of strongly correlated electron systems 
using this method. 

In the homogeneous case, the effectiveness of the method was demonstrated in Ref . [|J . A 
closely related approach also appears in Ref. where the gap equation in a superconductor 
has been derived using effective action ideas. In field theory this method also received some 
attention recently in treating a scalar field with a quartic interaction 0. Here we carry out 
the calculations of the energy in the non-homogeneous case as far as possible analytically. 
It will be seen that the method provides a very compact expression for the energy that goes 
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beyond the RPA method in a natural way. We also apply our result for the energy to the 
homogeneous case. We explicitly calculate the effect of taking exchange into account in 
the ring diagrams. Our results agree well with results obtained using quantum Monte Carlo 
methods. 

In the following, we will present a way to compute the free energy of the electron system. 
Zero temperature results follow easily from this calculation. The system is assumed to be 
in an external static potential V(x). The result that we get is very general and could be 
applied to many different systems. 

The paper is laid out as follows: In the next section, section 2, we introduce the thermo- 
dynamic potential Q for a system subject to external sources. Q therefore will be a functional 
of these sources. The functional Q is the generator of connected Green's functions at finite 
temperature. Usually, one point source functions are used. Here, we will instead introduce 
two-point functions @, ||, [pfl , fllQl . The introduction of two-point external sources enables 
us to take into account the higher order corrections to the stationary phase approximation 
of the partition function in a better way. By including merely two diagrams (expressed in 
terms of the variables which are conjugate to the external sources) in the effective action, 
we are able to obtain the contribution of all the second order diagrams (and beyond) of 
the Coulomb interaction to the energy in a very compact way (see Fig. [3]). This is the 
crucial advantage of this version of the effective action. The diagrams which are included 
here are two-particle irreducible, in contrast to the one-particle irreducible diagrams that we 
get when only one-point functions are used in the effective action. However, before doing 
this, we introduce a new field to replace the quartic Coulomb interaction by the well known 
Stratanovich-Hubbard transformation, so that all the interactions become local. This new 



field is simply the Hartree potential as in the RPA method [11]. However the treatment 
presented here is somewhat different from the one in the aforementioned reference. Here, we 
explicitly introduce a term that describes correlation of the Hartree field at two locations. 
The sources can also be taken to be instantaneous so we can get a time-translation invariant 
solution. 

In section 3, we introduce the finite temperature effective action T. We get this action by 
a Legendre transformation from Q. Therefore, we introduce new variables conjugate to the 
sources. This is exactly the same type of transformation as is used in thermodynamics. At 
the tree level, the T functional is simply the action of the system and is related to the free 
energy of the system. For homogeneous systems, T is known as the effective potential. In 
our case T will be a functional of three variables, the Hartree potential tp(f), the correlation 
function C(x,y) of the Hartree potential, and the Green's function p a/ 3(x,y) of the Fermi 
field. If the external sources are taken to be instantaneous, the points x = (ri,f{) and 
V = (T2,r^), correspond to the same temperature, T\ = r 2 , and the function p a /3(x,y) simply 
becomes the density matrix of the system, a and j3 are spin indices. Even though it is possible 
to express the free energy in terms of the density matrix, it makes the manipulations more 
cumbersome. 
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Figure 1: Graphs a, b, c and d are part of the expansion terms in T. Graph e is however 
reducible and does not show up in the expansion. The solid line represents the propagator 
p(x,y), the dashed line represents the propagator C(x,y). 
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In section 4, we solve for ip(x) and C(x, y) in terms of p a p(x, y). We are able to do this 
simply because in reality p a /3(x,y) is the only independent variable. In this case we get 
an expression for the effective action V in terms of the density only Therefore getting an 
expression for T in terms of p a fs{x,y) is the essential result of this work from which other 
important calculations can be made. 

In section 5, we give an explicit expression for the free energy in the homogeneous case 
in a neutral background. Here we go beyond the RPA results. We show that new diagrams 
appear at finite temperature. The famous 'anomalous' diagram is just one of them (see 
Fig. [|). An important outcome of this calculation is the disappearance of the InT behavior 
in the specific heat that is due to exchange |T3[ . 

In section 6, we treat the zero temperature case. We show explicitly that the anomalous 
diagrams disappear as they should. Here again we go beyond RPA to include exchange 
effects on the ring diagrams. We also show that some of the ladder terms are included in 
our approximation. 

In section 7, we continue the treatment started in the previous section and calculate the 
contribution of the ring diagrams if they include exchange. 

In the conclusion, we reexamine the non- homogeneous problem. We show how this 
method is related to the DFT formalism. In fact it is shown that a statement like the 
Hohenberg-Kohn theorem is trivially realized in this formalism. Similarly the question of 
v-representability can be given an affirmative answer within perturbation theory. 
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2 The Thermodynamic Potential Q 



In this section, we obtain an expression for the functional Q to 'order' h 2 . This functional is 
the logarithm of the partition functional Z of the system in the presence of external sources. 
The Z functional is written as an integral over all possible allowable paths of the different 
fields with a weighting factor that depends on the value of the action of the given path. 
The functional Q is simply the generating functional of connected Green's functions at finite 
temperature of the fields involved. However because the Coulomb interaction, quartic in the 
Fermion field, is hard to integrate, we introduce an auxiliary field by means of the well known 
Hubbard-Stratanovich transformation so that we can avoid the quartic interaction and its 
nonlocal behavior. In other words, the procedure consists of transforming the problem into 
one in which the electrons interact locally with a Hartree-type potential. In this treatment, 
exchange energy terms will show up directly because of the anti-commutation of the Fermi 
fields, which is built into the calculation as will be described below. The correlation terms 
are due to the quantum fluctuations around the Hartree potential. 

In the following we give an explicit derivation of the terms to order h 2 in Q and no 
counting will be necessary. What we mean here is that we are including effects of second 
order in the full Hartree field of the theory. The equation of motion satisfied by p(x, y) 
goes therefore beyond the Hartree-Fock approximation and involves second order exchange 
effects. 

A non-relativistic interacting electron gas in an external potential V(x) has the following 
Hamiltonian in the second quantized form: 



H n = / d 3 x 



{x)VH a {x) + V(x)¥ a (x)V a (x) 



■/■ 

+ \f d^y-^^Ux^Uy^M^aix). (1) 

^ j \ x y\ 



Here, we use units such that h = m = e = kB = l and j3 = ^. ^(x) is a two-component 
temperature-dependent electron field, a and u are spin indices, i.e., a = 1 for spin up and 
a = — 1 for spin down. Summation is implicit for repeated indices. In the following x will 
always mean a 3-D space vector. The system is constrained by the condition 



J d 3 x^ a (x)V a (x) = N (2) 

where N is the electron number operator which is constant. Therefore the density operator is 
simply ^ a (x)^ a (x). Because of this constraint, we prefer to work instead with the following 
Hamiltonian 
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H = H -nN 



(3) 



where \i is a Lagrangian multiplier. The electron operator can be expanded in terms 
of a complete orthonormal set of one-particle functions <fk(x), so we write 



*t( x ) = ^al^(x). 



(4) 



afe and a£ are annihilation and creation operators. The subscript k includes momentum and 
spin. The ground state and the excited states can be represented as Slater determinants 
formed by the wave functions <pk{x)- In our case, it will be advantageous to take these 
functions to be self-consistent Hartree eigenfunctions. In the homogeneous case they reduce 
to plane waves. Since we are going to use a path integral formulation, we give the Lagrangian 
associated with the Hamiltonian H. 



L = J (f'.r 
1 



d 1 

M\x)—ty(x) + -ty(x)V 2 ^ j (x) - V(x)^(x)^(x) + /2^(x)^f(x) 



. d 3 xd 3 y &(x)&{y)V(y)V{x). 
* .1 \x y\ 



(5) 



So for a system at thermal equilibrium, we have to replace the time component t in the 
Hamiltonian by — it. The partition function Z is as usual defined as a sum over all possible 
states <£>. In the absence of external sources, we have 



Z[0] = ^<$|e^|$> 

= nJ DVD&exp (S[V,&]) (6) 

where N is a normalization constant and the integration measure of the Fermi fields is only 
defined for fields that satisfy: 

y(0,x) = -#(/?, a). (7) 

S[$, v&t] is the action of the system 
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drL 



= J dr J d 3 x ^(x 

I ft 3 r i 
- - / dr d 3 xd 3 y — — z;W(x)W(y)*(y)*(x) 

* Jo J \ x ~ y\ 



(8) 



The functional integral for the partition function becomes 



Z[0] 



N J DtfDtf+expj- J dr x d 3 x &(x) (J^ - ^ V 2 - n + V(x)^j 

4 /V*-, f^'l^y Sl r-T,) 
2 Jo 4/ F-2/1 



(9) 



In writing the last term we made use of the fact that the Fermi fields satisfy a Grassmann 
algebra. An infinite self-energy term has been dropped from the above expression. Such a 
term cancels at the end and has no effect on the evaluation of the eigenfunctions or related 
physical quantities. Now, we set for notational convenience: 



G-\x, y) = (£r ~ " f + s ( x ~ v) 



(10) 



and 



A(x,y) = 



x r y 



x-y 



(11) 



Before proceeding further, we have to get rid of the quartic term in the Lagrangian. This is 
done by introducing an auxiliary boson field. 
We write 



Z[0] =N' J DVD^Dip exp d 4 xd 4 y *t (x)G _1 (a;, y)* Q (y) 

- l - J d*xd 4 y ip(x)A~ 1 (x, y)ip(y) + J d 4 x <p(x)*t,(x)*(x)} . (12) 
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The fourth component is the time component integrated over the proper range. It is easy to 
see that by using the following formula, 



J dx e -^ x+bx = (det Q)"5 e^ Q ' lb (13) 

and integrating over (p we get back the original expression for Z[0]. The prefactor N' is a 
new normalization constant. The operator A(x, y) is clearly invertible, 

A- 1 (x,y) = ~S(x-y)V 2 . (14) 

The non-local character of the Coulomb interaction has been removed by introducing the 
new bosonic field tp(x). It can be shown that <p(x) is the Hartree potential by using the new 
equations of motion derived from the new action of the problem. We also couple the fields 
to local and non-local sources J(x),Q(x,y) and B(x,y). Now we introduce the functional 
Q, a generator of connected Green functions. This is introduced through the normalized 
partition functional which is now a functional of the external sources. It is given by: 



Z[J,B,Q] = exp{-pn[J,B,Q]} (15) 



such that 



exp(-Q[J,B,Q]) J DVD&Dip exp (-£[#, * f , <p]) = 

J D^D^D^exp!^-S[^,^,^] + J d 4 xd 4 y^(x)Q( y x,y)^(y) 

+ d 4 xd 4 y<p(x)B(x,y)ip{y) + J d 4 x<p(x)J(x)} (16) 

where we have set (3 = 1. S[^f, ^\ <p] is simply the action of the transformed problem, 

= j d 4 xd 4 y ^(x)G-\x,y)^(y) + ^(x)A-\x,yMy) 

- J d 4 x(p(x)^(x)V(x) (17) 
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Now we define three new variables: 



6J ( x j \J=B=Q=0 = {<p(x)) = (p c (x) (18) 



S ~M(x 1 ^y \ J=B=Q=0 = \ ^Mv^ = \ [^(^)Vc(y) + C(x, y)\ (19) 
Sn[J,B,Q] 



6Q(x,y) 



j=b=q=o = (*Ux)*{)(v)) = Pap(x, y) (20) 



ip c (x) is the expectation value of the field tp(x) in the ground state. C(x, y) is a correlation 
function of the field <p{x). p a p(x,y) is the Green function of the Fermi field. p a p(x,x) is 
therefore the density of the system. Here and below, the "time" ordering operator is not 
written explicitly. Therefore terms like p{x, y) are defined by setting r x — r y = + . Note that 
C(x, y) measures the departure from quasi-independence due to the correlation between the 
values of the potential at two different locations. The expectation value of the Fermi field 
ty(x) is zero. In the following, we will obtain an explicit expression for Q[J, B, Q] to 'order' 
h 2 . Then by solving for (p c (x),C(x,y) and p(x,y) in terms of J(x), B(x,y) and Q(x,y) we 
get an expression for the effective action T which is a functional of the new variables. This 
will be done in the next section. 

Now we turn to determining Q. First we expand the exponent in the above integral around 
— \[rt — o anc i ip — ip Q where <p is the configuration of (p that extremizes the action S. 
Therefore, we have 

(A- 1 + B)i Po = -J. (21) 

It is understood from the above that there is an integration over space and time on the 
L.H.S. of this expression. We choose from now on not to write integrals explicitly unless 
there might be some confusion. Now we expand around (po, so we write 



P(old) = P(new) + V?0- (22) 



Then assuming that the main contribution to the integral comes from the saddle point, we 
get the following expression for the partition functional: 
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x exp 



detfC" 1 ] det^- 1 ]-^ 



where we have defined S and E to be 



1 



&(G- 1 + Q)*--<p(A- 1 + B)<p 



a = exp 

£ = J DVD&Dip exp 



1 



In the above, we have used the fact that 



da da) e~ atMa = det M 
for Grassmann numbers a and a*, and 

J dx e ~^ xMx = (detM)- 1 * 
for c-numbers x. The argument of the first term on the right is 



s [<Po] = \p* A Vo + ^PoBtpo + Jtfio 

and 



G' 1 = G- 1 - Vo - 



Now using the fact that 



detA = e TrhlA , 
we can get an explicit expression for Q[J, B, Q\. 
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Q[J, B, Q] = S(ip ) - Tr ln[l + G y? + GQ] + \ Tr \n[l + A B] 

-- [dVD&D<p ((ptfV) 2 ^ + ... (31) 
2 J L 

After finding Q, we now calculate the effective action T. We have to solve for the external 
sources in terms of the physical variables p, C and y? c . This we do in the next section where 
we find the effective action at finite temperature. However, this treatment applies to zero 
temperatures too. In this case too, we get an expansion in h . Again h is not a true expansion 
parameter in the true sense of the word. It is only used as a bookkeeping method for the 
diagrams included in the expansion of Q. 
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3 The Effective Action T 



The effective action is a functional of ip c (%), C(x,y), and p(x,y), which is obtained by a 
triple Legendre transformation from fl[J, B, Q] 

V[(p c ,C,p] = Q[J,B,Q] - J d 4 x <p c (x)J(x) - ^ y d [xy] ip c (x)B(x,y)ip c (y) 

J d[xy] C(x, y)B(x, y) - J d [xy]p(x, y)Q(x, y) . (32) 



1 

2 



where we have written d [xyz..] for d A xd 4 yd A z... to simplify the notation. It is understood 
that the integration over time is carried for r in [0, H0\ for the nonzero temperature case 
and over all time for the zero temperature case. In this section we will assume we are in 
the zero temperature regime. However this discussion applies equally well to the nonzero 
temperature case. 

It is easily verified that : 



6T[ip c ,C,p] 
8<p c (x) 



= -J(x)- J d 4 y B(x,y)<p c (y) (33) 



5C{x,y) 2 

and 

5T[ Vc ,C,p] 



ST[ ^ p] --^y), (34) 



-Q(x, y) . (35) 
8p(x,y) 



When we turn off the external sources, the above equations give the values of (p c ,C and p 
that minimize the effective action. <p c (x) and C(x,y) are really dependent variables since 
they depend on the auxiliary field <p(x), so we should be able in principle to express them 
in terms of p(x, y) which gives us the density. 

Now we have to express J, B and Q in terms of ip c , C and p. First we note that when 
h — 0, we have 



Tj = ^ (36) 



Hence, we can write 
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Po = <p c + P (37) 
where tp is of order h. Similarly, we write 

S(<p )=S(<p c )+hS 1 . (38) 
Now we seek an expression for T in the form 

r = r + nr 1 + ^ 2 r 2 . (39) 

It then follows that 

r = ^ C A-Vc. (40) 

Now, we find approximate expressions for B and Q in terms of <p c , C and p. Using Eq.(31) 
and Eq.(35), we get the following relation: 

p~ x G = 1 + GQ- eGp c + 0(h) . (41) 

The expression for B is more involved. Again using Eq.(35) and treating ip as a functional 
of J and B, we get: 

BC = (l-2epp c ) (l-2^j -A- X C. (42) 
Inserting back all these expressions in T, we get for I\ 



r\ = -Tr In (p^G) + ^Tr ln(C _:L A) + ^Tr {A~ l C) - Tr {pG' 1 ) . (43) 

The terms of order h 2 are more involved to treat but again straightforward. The steps are 
similar to those in the case of the effective action with one-point sources only [[L2| . Actually, 
it can be shown that the next terms in r 2 are the sum of two-particle irreducible diagrams 
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of the theory ||,||. Here, we sum only the first two diagrams in this series expansion. 
Fig. m shows the first four diagrams in this expansion. The fifth diagram is not part of the 
asymptotic expansion of T. These first two diagrams enable us to include first and second 
order exchange effects. Diagram e is not part of r 2 since it is two-particle reducible. However 
it is one-particle irreducible and it is part of the usual effective action fill , 



r 2 = -^e 2 J d[xy] p(x,y)p(y,x)C(x,y) 

- ^e 4 J d[xyuv] p(x, y)p(y, u)p(u, v)p(v, x)C(x, u)C(y, v). (44) 

Hence the reducible graphs do not appear in the expansion. It is the term S that appeared 
in Eq.(23) that is responsible for excluding such a graph from the expansion. T 2 and higher 
order terms represent vacuum graphs where the propagators are p(x,y) and C(x,y) of the 
Fermi field and the boson field respectively. Therefore to order h 2 , the effective action has 
the following expression 

r[<p e ,C,p] = -S(<p e ) + ~Tr InAC-' + A-'C (45) 
-Tr \np- 1 G + G- 1 (i Pc )p + T 2 [C,p}, 



where 



G-\ Vc ) =(j^-\^ 2 -^ + V{x) - <p e {x)) . (46) 

The term p J d 3 x p(x, x) cancels the term pN in the Hamiltonian H, and therefore it can be 
handled without difficulty. However, we have to deal with the p that appears in the term 
Tr In G. This could be removed by a convenient choice of the path of integration in the 
complex a;-plane. The case p = occurs when there are no charges present. We are dealing 
with negatively charged electrons bound by a static potential V(x). Therefore bound states 
appear, with negative energies bounded from below for any physically sensible V(x). In 
Fig. H we show how by going from the C\ path, which corresponds to nonzero p, to the Co 
path we pick up contributions from the bound states between the two paths, since energies 
of bound states are poles of the propagator G^ 1 in the complex u plane. Note that the 
eigenvalues u>k of these states do not take into account correlations due the Coulomb field. 
Because of the logarithmic operator, there is a cut along the real positive axis. By going 
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back to real time and then Fourier transforming the time dependence, we get 
J dt Tr In - i V 2 - p + V(x)^j 

=l dt L^i tri «{-"-\ v * +v ) + ^ (47) 

where tr applies only to spatial variables. The u k s are the eigenvalues of the following 
equation 



-iv 2 + V{x)j <p k (x) = u k ^ k {x) . (48) 

The single particle functions <p k (x) of the potential V(x) are not supposed to be used alone 
as the starting point for any numerical calculations since they do not take account of the 
repulsion between electrons. Any numerical calculations will have to start by solving Eq.(56). 
However, it is clearly advantageous to separate the effect of the external potential V(x). This 
separation also appears in the usual treatments. The ground state energy E g is given in terms 
of the effective action per unit time. Using the above results, and the fact that the time 
independent values of p(x,y),C(x,y) and ip c (x) are defined by setting t x = r y = 0, the 
expression for E g is 

E g J dt = ^ J d[xy] Vc {x)A-\x,y)ip c {y) + ^Tr \nC A' 1 
-\\ d[xy] A' l (x,y)C(y,x) - ^ u k j dt 

-Tr XnpG-'-r j d'x ( A _ I V 2 + V(x) - e <p c (x)j p(x,y)\ x =l Ty 
+ ^e 4 J d[xyuv] p(x, y)p(y, u)p{u, v)p(v, x)C(x, u)C(y, v) 

+ ^e 2 J d[xy] C(x, y)p{x, y)p(y, x) (49) 

where we have dropped terms of order higher than e 4 in Y2 for simplicity. This is the full 
expression for the ground state energy of the system. Here G^ 1 has p set to zero. To be able 
to use this expression for E g , more approximations are required. In later sections, we use 
this expression to find the energy of a homogeneous electron gas both at zero and nonzero 
temperature. Since tp c is a dependent field, then in principle we should be able to solve for (p c 
and C in terms of p only. This we do next but only approximately to keep the calculations 
manageable. 
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Figure 2: Path of integration used to obtain Eq.(47). 
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4 The Effective Action as a Functional of p(x, y) 



In this section, we find an expression for T solely in terms of the Green's function p(x,y). 
This is done by finding the expressions for ip c (x) and C(x,y), that minimize T, in terms of 
the density p(x,y). For (p c (x) we get the following expression: 



<Pc(x) = -e J d 3 y j| 



p(y,y)_ 
y\ 



Similarly, minimizing T with respect to p(x,y), we get 



(50) 



S(x, z) = 



JL. _ i V 2 + V(x) - e <p c (x) 



e 4 / d[yuv]p(v,u)p(u,x)p(y,v)p(z,y)C(u,y)C(v,x) . 



(52) 



Finally, minimizing T with respect to C, we get 



6(x,z)=J d 4 y[A \x,y) - e 2 p(x,y)p(y,x)]C(y,z) 



-e 4 / d[yvu]C(u,v)C(y,z)p(x,u)p(u,y)p(y,v)p(v,x) . 



(53) 



To be able to solve for C(x,y) we have to linearize the theory. In the following, we keep 
terms only to 'order' e 4 . Hence, we find that the correlation C(x, y) is given by 



C(x, z) = A(x, z) + e 2 / <i[xix 2 ] A(x, X\)D{xi ) x 2 )A(x 2 , z) 



+e 4 J d[x 1 x 2 x' 1 x' 2 ] p(xi, x' 1 )p(x' 1 , x 2 )p(x 2 , x' 2 )p(x' 2 , x 1 )A(x 2 , x 1 )A(x' 2 , z)A(x, x[) (54) 
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Figure 3: The expansion of the propagator C(x, z) to order e 4 . 



where we have set 

D(x 1 ,x 2 ) = p(xi,x 2 )p(x 2 ,xi) . (55) 

The above equation for C(x, z) can be represented diagrammatically as shown in Fig. [3[ 
The first term is the bare Coulomb potential. The second and the third are the direct and 
exchange term with one bubble. So far the Fermion propagators are the true propagators 
in this expansion. In the following we will ignore the corrections with two bubbles and 
higher. Clearly C(x, y) is the full screened Coulomb potential if all orders of e 2 are included. 
Similarly, to order e 4 , the equation satisfied by p(x, z) is 



^~2 V +V{x) + e J dy W^\ 



p(x, z) = 5(x, z) 



-e 2 / dyA(x,y)p(x,y)p(y,z)-e 4 / d[yxix 2 ] p(x,y)p(y, z)A(x,x 1 )D(x 1 ,x 2 )A(x 2 ,y) 



+e 4 / d[yxix 2 ] p(x 2 ,xi)p(xi,y)p(x,x 2 )p(y, z)A(x 1 ,x)A(x 2 ,y) . 



(56) 
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Figure 4: Born approximation for two-particle Green's function. 



This is obviously the known equation satisfied by the Fermi Green's functions [11]. The 



term on the R.H.S. of order e 2 is an exchange term or Fock term due to the statistics of the 
electrons. The last two terms of order e 4 take into account collisions and are equivalent to the 
usual Born approximation for two-particle Green's function in scattering theory, Fig. [|. The 
most likely way to solve these integro-differential equations is by iteration. Using the above 
equations, we can write an explicit expression for the energy in terms of the full propagator 
p(x, z) of the theory. However, such an expression suffers from the problem of over-counting 
some of the states of the system. This is mainly due to the fact that e 2 is not a good 
expansion parameter. Hence physical arguments are used to discard some terms at this level 
of approximation and instead include them at higher orders. An expression for the electron 
propagator in terms of the free propagator is given in Fig. [5] . The corresponding analytical 
expression can be easily written following usual rules [TT| . Here we keep only the first three 



terms. The second term vanishes in the homogeneous case, the case we are interested in this 
paper. We expect the expression for the energy, which has not been linearized, to be a good 
one if we believe that a stationary phase approximation is viable. Without including T 2 , the 
expression obtained for the energy is correct with p(x, y) the Hartree propagator and (p c (x) 
the Hartree potential. This indicates that a stationary phase approximation is possible and 
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Figure 5: Approximate solution to the one-particle Green's function. The propagators on 
the right correspond to the free theory. 
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that the higher order correction T 2 should be a small perturbation to the Hartree solution. 
The linearization of the problem must take account of this. Hence, we can justify the validity 
of the expansion of T in h, 



r = T + hT 1 + h 2 T 2 + ... (57) 

by the result we get in the end. To second order in h, Eq.(54) becomes 

J dyA~\x, y) C(y, z) = S(x, z) + he 2 J dyD(x, y)A(y, z) 

+he A I d[yuv] p(x,u)p(u,y)p(y,v)p(v,x)A(u,v)A(y, z) 



where p(x, y) is the free electron propagator. We have used only a first order approximation 
to the true propagator to find this equation for C(x,y). Using this equation, we get an 
expression for the energy to order h 2 with proper symmetry factors for the diagrams involved 
in the expansion. Setting h — 1, we obtain for an interval of time T, 

T EM = 6 - [ dxdy P"( x >*)P»(y>y) + e f dx V(x)p H (x, x) 
2 J \x-y\ J 



+ 



/ dx ( -^V 2 p H {x,y) J \ x=y -T^2u k 
J V 1 ' ^ 



\ g4 / dxi dx ^ dxd V PH(xi,x)p H (x,x 2 )pH(x2,y)pH(y,x 1 )A(x,y)A(x 2 ,x 1 ) 



H — Tr In 



e I dydx x dx 2 p H {x u x)p H (x, x 2 )p H (x 2 , y)p H (y, x l )A(x 2 , x 1 )A{y, z) 
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2e 4 / dydx 1 dx 2 PH(x,x 1 )p H {xi,x 2 )pH(x2,y)pH(y,x)A(x2,xi)A(y,z) 



-Tr In 



S(x,z)-e 2 J dy p H (y,y)p H ( y x, z)A(x,y) - e 2 J dy A(x,y)p H ( y x,y)p H ( y y, z) 
-e 4 J dydxidx 2 Ph(x, y)pn(y, z)p H (xi, x 2 )ph(x 2 , x 1 )A(x 1 ,x)A(x 2 , y) 



+e 4 J dydxidx 2 p H (x 2 , x 1 )p H (x 1 ,y)p H (x, x 2 )p H (y, z)A(x 1 ,x)A(x 2 , y) 



-e 2 dy p H (y,y)p H (x,z)A(y,x)\ x=z + e 2 / dxdyA(x,y)p H (x,y)p H (y,z) 



-e 4 J dxdydxidx 2 8(x - z) A(x, xi)A(y, x 2 )p H (y, z)p H (x 2 , x 1 )p H (x 1 ,y)p H (x, x 2 ) 

+e 4 J dxdydxidx 2 S(x - z)p H (x, x 1 )p H (x 1 , x 2 )ph(x 2 , v)ph{v-, x)A(x 2 , xi)A(y, z). (59) 

The trace, Tr, acts on x and z. So now we have obtained an expression for the energy 
in terms of the Green's functions p(x, z) only. This is a major result of the work. As we 
mentioned above, we solve for p(x, y) to order e 4 by iteration of Eqs. (52-53). The expression 
for the energy can be solely written in terms of p(x, x). However as we stated earlier this is 
not advantageous and it is better to keep working in terms of p(x,y). 

We see that the classical Coulomb term and the interaction with the external potential 
appear naturally in this approximation. They can be easily separated from the full expression 
for the energy. Another important point is that a gradient of the density also appears 
naturally within the above expressions. Therefore the kinetic term is treated more accurately 
in this method compared to the density functional method. By expanding the first logarithm, 
we immediately obtain the Hartree-Fock exchange term, i.e., 



t? 1 2 / . , pH(x,y)p H (y,x) 

E exch = -e / dxdy — — 5(t x - r„) 

2 J \x-y\ 



(60) 



The usual expression for the exchange energy follows trivially from the previous expression 
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upon using the commutation relation of the Fermi field. The infinity arising from this cancels 
the self-interaction term that was dropped in Eq.(9). There is a similar exchange term that 
comes from the second logarithmic term . This term is canceled by another similar term 
outside the logarithm. The usual RPA term is contained in the first logarithm. The new, 
extra term in the same logarithm provides among other things exchange corrections to the 
ring diagrams. The higher order exchange diagrams are also included in this approximation. 
We will say more on this when we treat the homogeneous case in the following sections. 
Before we go on to the next section, we should point out that an expression for the energy 



in terms of the density only is possible in this formalism [ 13| . The expression is complicated 



however and does not seem to be useful. We were also able to show that a derivative 
expansion in terms of the Hartree field is more suitable in our case rather than seeking one 
in terms of the density. In fact a derivative expansion in terms of the Hartree field seems 
more appropriate since it is better behaved than the density for problems where the density 
changes sharply. We hope to address these questions in a later communication. 
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5 The Fermi Gas at Low Temperature 



In this section we use the result found previously to improve on calculations of the exchange 
energy and ring diagrams made by Isihara and coworkers and others and calculate the specific 



heat at low temperature |L4| . A short review of these calculations appears in Ref.|15| and 
offers a good comparison of all these methods up to the year 1975. The authors of this 
paper end with an interesting conjecture which we prove here to be true. The appearance 
of temperature dependent logarithmic terms in the internal energy found through various 
calculations is unacceptable. These terms will reappear again in the specific heat and will 
spoil the observed linear behavior at low temperature. In the past these difficulties were 
not given close attention. They were either hidden in an effective mass or the expansion 
was not completely correct as we argue below. The history of this problem dates back to 
the '30s when it was first pointed out by Bardeen in Ref. [T5| . A preliminary treatment of 



this work appeared in Ref. Q]. It was argued in Ref.[(|] that a correct approximation to the 
internal energy of a many-particle system up to exchange, i.e., a Hartree-Fock approximation, 
will automatically include summing an infinite set of diagrams generated by the known 
"anomalous" diagram, Fig. M. 



At low temperature, the ring diagrams treated in Ref.|]l7| are of order e 4 . The exchange 
energy is of order e 2 . However at nonzero temperature, we have another contribution to 
order e 2 that comes from summing an infinite set of diagrams that were never treated before. 
Therefore the new set can not be ignored in any finite temperature treatment. We show below 
that taking these diagrams into account is exactly what is needed to solve the anomaly in 
the exchange energy at finite temperature. The grand partition function Z of the system is 
given in this case by the following expression: 



Z = exp(— (3Vt) = Trexp I —(5 



[R - fiN) + j<p + + ^Qi/j 



(61) 



The thermodynamical potential is a functional of j, B and Q but it is a function of the 
variables p, T and V. To get the free energy F[p, <p c , A] we make the following triple 
Legendre transformation 



F[p, <p c , A] - p N = B, Q] - jip c - -ip c B<p c -pQ-AB 



(62) 



where p, tp c and A are defined as before, Eqs. (18-20). We have seen that F (E at T = 0) 
can be calculated perturbatively by expanding the exponent around the classical Hartree 
potential in a neutral background, i.e., with no external potential. The coefficients of this 
expansion are expressible in terms of Feynman diagrams. Here in this section we keep only 
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Figure 6: A new infinite set of diagrams at nonzero temperature. 
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the first two terms of the expansion, in other words we only include the first diagram in 
Fig- 0- 

The equation of motion satisfied by p(x, x') is given by: 



' d_ _ 1 

lh~2 



where A(z — x) is the bare Coulomb potential to first approximation. It must be pointed out 
again that a full solution to A(z — x) gives the expected shielded potential. An approximate 



solution to the full linearized equation for A(z — x) was first given in Ref. ||18||. We must 
stress however that in what follows we use only the zero order solution to A(z — x). In this 
case the above equation is the Hartree-Fock approximation to the equation of motion for the 
one-particle Green's function p(x — y) at finite temperature |[L9| . The final expression that 
we are after is that for the thermodynamic potential Q(p,T, V). Hence within the above 
approximation, the expression for the thermodynamic potential has the following form: 

Q = n + n ex + ^Tr {In [S(x - y) + e 2 ^(x - y)] - e 2 ^{x - y)} 

- Tr {In [8{x - y) - e 2 5s(x - y)] + e 2 5s(x - y) } (64) 

where the functions 3?(x — y) and 5s(x — y) are given by 

P 

'3, 







dt(x — y) = J dr J d zA(y — z)p(x — y)p(z — x) (65) 

and 



rP 

5s (x — y) = I dr I d 6 zA(x — z)p(z — x)p{y — z) . (66) 



o 



The term Q is the free contribution, Q ex is the usual first order exchange term, the third 
term represents the usual ring diagrams and the last term represents a new series of diagrams 
shown in Fig. The function A(x — y) is the bare Coulomb potential and p(x — y) is the 
free one-particle Green's function at finite temperature. It should be stressed from the above 
treatment that the new diagrams are necessary for the consistency of our approximations. 
The first diagram of the new set of diagrams was treated previously by Isihara and Kojima 



in Ref.pOj in a calculation of the specific heat. Here we do the same calculation again 
but taking into account the full sequence of the new diagrams. The method used by these 
authors differs from ours, however they are both based on path integrals and both start 
from a partition function. In this sense our methods are more exact than those that start 
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from a zero temperature treatment like Gell-Mann's treatment in Ref.pTfl or those based on 



a phenomenological theory like Landau's. The authors of Ref. |TJ] claim at the end of their 
paper that screening is probably not needed to have a reasonable answer for the specific 
heat. Here we show that indeed this claim is true. It is well known that the exchange energy 
has an odd InT behavior that ends up contributing a TlnT term to the specific heat. This 
is obviously not reasonable for low temperatures. The exchange energy at low temperature 
is given by [FM : 



a 



-v 



e 2 m 2 
2tt 2 



7T 



-a 



TV 



3/3 2 



l+ln/fy-7) 



(67) 



Now we give the contribution due to the ring diagrams at finite temperature ||23|| . This 
contribution is given by 



1 



R = -Tr {In [S(x - y) + e 2 $t(x - y)] - e 2 $l(x - y)} 



(68) 



In momentum representation, this is given by 



y n ~°? Air r°° 

W n fr^ ( 27 3 Jo 



(69) 



where p represents momentum and u n = n = 0, ±1, ±2, .... The function 5R(p, oj n ) 

is given by: 



dq- 



W+q) 2 -W 



WP 2 J - ^l+^ip + q) 2 -^) 
The function f(q) is the Fermi-Dirac function, 



-(f(q)-fip + q))- 



(70) 



1 + exp 



P( q 4 -A*) 



(71) 



Therefore the finite temperature contribution of the ring diagrams to leading order in 
r] = (3fi, and for high density is given by 



R = AT(0.0622 lnr s - 0.142) + 0.0181 N(T) \i 



,ar s {T), 



(72) 
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Here N is the total number of particles, which is a constant. All the other variables are 
defined as follows: 



/ 4 ,1 me 2 h 2 „ T/rTn . V ,„ . 

° = <S>*- r ' (r) = 5m? = S?' W(T) = f^rR • < 73 > 

To obtain this answer, we have used explicitly the fact that r s is small to simplify the 
calculation. It is however important to stress that this approximation does not affect our 
final result regarding the InT terms in the specific heat. Obviously our result will not apply 
even to r s ~ 1. 

The contribution of the new set of diagrams, Fig. 0, is given by : 



A = Tr | In [ 8(x — y) — e 2 Q(x — y)] + e 2 Q(x — y)\ 

4E -T% / x " dx I ln i 1 - e ^ "»)) + e ^ ^)] ■ (74) 

P n=-oo ( 27T ) JO 



■r- is the normalize^ »..„,>„.,+ »>,, • — '-" 1 - |J " 

function $S(x,u n ) is given by : 



Here x = is the normalized momentum and u n = - ^ for n = 0,±1,±2, .... The 



%(x,u n ) = ^p(x,u n ) I d 3 p- -^— 2 — — ^ (75) 

1 + exp(?7(^- — IJJ | p — xkf I 



7T 2 ' 



with 



p(k,u n ) = - — 2 (76) 

iu n - (— - (l) 

and r\ = P/j,. The mass of the electron is taken to be m=l. The calculation is carried out for 
large rj, i.e. for low temperature. Again here we used explicitly that r s is a small parameter. 
After some lengthy calculations, we get the following explicit expression for A to order -y : 

where £ is a constant. 
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There are two things to note. First there is no zero temperature contribution due to these 
diagrams, which is as it is supposed to be. Second, the sum of the new diagrams gives an 
answer which is of order e 2 , that is it has the same order as the first order exchange energy 
even though the terms in the sum are of order at least e 4 . This is the reason we believe that 
an expansion in e 2 did not pick up this new infinite set of diagrams. It is also of interest to 
note that the ring diagrams add up to a term of order e 4 . Now a close comparison of the 
three terms shows that the In T term in the exchange energy gets cancelled by the one found 
from the new set of diagrams. The ring diagrams contribute a temperature dependent term 
of order e 4 and hence they are not needed to solve the InT-term in the exchange. In other 
words screening does not play any role in this odd behavior of the exchange energy or for 
that matter of the specific heat as we show next. 

The specific heat formula is: 



-T 



,d 2 F 
gf2 



(7E 



where the free energy F is given by, F(T, V, N) = Q + /iiV. Using standard techniques p4|1 . 
we calculate the contribution of all of the above terms to the specific heat. To first order in 
■h, we have |23[ : 



k 2 T 2 

i4 



F(T,V,N) = F (T,V,N) + 



20.19 



0.077 lnr s + 0.184 



(79) 



where Fq is the free energy of a free electron gas. Therefore, if C v is the specific heat for 
the ideal Fermi gas, we get for the specific heat at low temperatures and high densities the 
expression, 



= 1.00 + 0.185r s -r 2 (0.0017- 0.00071nr s ) . (80) 

The constants £ and 7 are approximately equal to 8.0 and 2.0, respectively. Now we turn 
again to the zero temperature case and calculate the correction due to the inclusion of 
diagram b, Fig. [I], in the expansion of the energy. 
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6 The Homogeneous Electron Gas at Zero 
Temperature 

In this section we apply the main result of section 4, the expression for the energy, to the 
homogeneous case at zero temperature. We assume that there is a background of positive 
charge of equal magnitude to the average density of the electron gas. Hence the system is 
neutral. The literature on this problem of the calculation of the correlation energy is huge. 
The most complete treatment, so far as we are aware, was given by Bishop and Luhrmann 



and it was restricted to zero temperature effects |25| . Since the system is homogeneous, the 
final expression for the energy will be given in momentum space. The Green's function that 
will be used in the following is the solution to the first iteration of the nonlinear equation 
satisfied by p(x,y). Since there is no external potential, the input Green's function is that 
of a free electron gas. The energy expression is given in imaginary time, hence we use the 
following expression for the free Green's function 



Po(w,k)= . i • (81) 

In this section and the next, the electron Green's function definition differs by a factor of 
i from the one given before. In terms of po(x,y), the ground state energy is given by the 
following, omitting the subscript for simplicity, 



-TE gs = J d*x QiVVfol/)) 



x=y 



e 4 f e 2 
+ — / dx dz 5{x — z) Ai(x, z) + — / dx dy p(x,y)p(y,x)A(x,y) 



~Tr^, g ) In [5{x,z) - e 2 H{x,z) + e A M{x,z) - 2 e 4 U(x,z) ] 



+Tr( x>z j In [8(x, z) — e 2 JC(x, z) — e 4 S(x, z) ] 

f fe 2 e 4 1 

+ / dxdzS(x-z) <— TC(x,z) — —Ai(x, z) + e 4 S(x, z) + e 4 II(x, z) > 



(82) 
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where 

M(x,z)= / dydx 1 dx2p(y,x 2 )p(x2,x)p(x,x 1 )p(xi,y)A(xi,x 2 )A(y,z) , (83) 



H(x,z) = J dyA(z,y)D(x,y) , (84) 



U{x,z) = J dydx 1 dx 2 p(y,x)p(x,x 1 )p(x 1 ,x 2 )p(x 2 ,y)A(x 1 ,x 2 )A(y,z) , 



(85) 



K,(x,z) = J dyA(x,y)p(x,y)p(y,z) , (86) 
S(x, z) — z) - S 2 (x, z) , (87) 
£i(x,z)= / dydxi<te 2 p(»,y)p(y,^)p(»i,»2)p(aJ2,aJi)>l(a!,iCi)>l(aJ2,y) , (88) 



Y, 2 {x,z)= J dydx 1 dx 2 p(x,x 2 )p(x 2 ,x 1 )p(x 1 ,y)p(y,z)A(x 1 ,x)A(x 2 ,y) . (89) 

The above equation, Eq.(82), is simply Eq.(59) taking into account of the fact that p(x,x) 
gets canceled by an equal and opposite charge. Because of the transitional invariance of 
the system, we can write a simple expression for the correlation energy in the momentum 
representation, 

Ec = j \\ , * 2 ^ A (Pl - P2)A(p 2 - P3)p(Pl)p(P2)p(P3)p(Pl ~P2+ Ps) 

1=1 

+ \v J -0^^ {In [1 + e 2 H(p) + e 4 M(p) - 2e 4 n(p)] - e 2 H(p) - e 4 M(p) + 2e 4 n(p)} 

" V J % { ln t 1 - ^(P) - e4E (P)l + + ^ } • ( 9 °) 

Since the ring diagrams are part of i? c , this will be the first test to see if our expansion gives 
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the correct leading result. First we show how this equation follows from Eq.(82). We start 
by letting 

H(x, z) = - J dy D(x, y)A(y, z) (91) 
or, since the system is homogeneous, we can write instead 



H(x — z) = J dy D(x — y) A(y — z) 



(92) 



Next, we rewrite D(x — y) and A(x — z) in terms of their corresponding Fourier transforms, 
i.e., 

/dk 
— exp[ik-(x-y)]D(k) (93) 

and 

My ~ *) = J J^y[ exp[ig • (y - z)]A(q) . (94) 
Hence, by convolution we have, 

H(p) = (A*D) (p) , (95) 

however since 

D{P) = J-^y 4 p(p + q) P (Q), (96) 

it follows then that, 

U{p) = A(p) J j^y 4 p( P + q)p(q). (97) 
After using the expressions for the propagators, we end up with the following expression for 

nip) 
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47r /" 1 1 

H[P} = fj (2^ -iu + ljp-p _ 2 ( tD + w ) + I(p + ? -)2_^ 

with q = (u,q). Integrating first over uj in the complex plane, we find that 

P 2 i(27T) 3 w 2 + (l( ?+5 )2_lg a )=» 

where we added a factor of 2 to account for spin. This term 7i(p) will prove to be all that 
is needed to reproduce the known RPA result. All other terms that appear in the energy 
expression are new additions to the correlation. Before showing this explicitly , we give the 
expressions for the remaining terms K,(p), A4(p), H(p), S 2 (p) and 

m=M J-^pip+^MQ) (ioo) 



M(p) = A(p) J * * p{kMh + p) P (k 2 )p(h - p)A(h -k 2 +p), (101) 
n(p) = A(p) | ^ JjLpthMh - p)p(h -p- q)A(q) , (102) 
£i(p) = p(p) J -^^-^zp(Pi)p(P2)p(p2 +Pi -p)A(p 2 -p)A{p-p 1 ) , (103) 

^4 A2 (Pi)p(P ~ Pi)p(Pi + P2)p(p 2 ) (104) 

and 

S(p) = / ^^^MW** + A* - p) - p) 2 - A(p - k 2 )A(p - h)} . 

(105) 

The Fourier transform A(p) is given as usual by 
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Figure 7: Second order exchange diagram in Eq. |90| . 
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Figure 8: Ring diagrams generated by H(p) in Eq. (pl|). 
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A{p\ 



47T 
I P | 2 



(106) 



The terms in Eq. (^0|) have the following meaning after we expand the ln-terms. The first 
term represents the second order exchange term, Fig. 0. The term TC(p) is responsible for 
generating the ring diagrams, Fig. [| The term Ai(p) is responsible for generating some of 
the ring diagrams with exchange. The cross terms give the remaining ring diagrams with 
exchange, Fig. |9]. The H(p) term generates ring diagrams with a self-energy insertion in each 
ring, fig. ITUl Hence the terms in the first In should contribute the most. Some of the terms 
that appear upon expansion of the second In term are shown in Fig. [11]. The term JC(p) 



generates the so called 'anomalous' diagrams, Fig. || The term generates ladder terms 



like those in Figs. [12] and 13. S 2 (p) generates terms like those in Fig. [14]. Cross terms of 



the last two terms are shown in Fig. [15]. Diagrams that appear in Fig. [16] are important in 
a nonhomogeneous medium and must be accounted for since they no longer vanish. They 
arise whenever we have a p(x,x) term in the expression for energy, Eq.(59). Normalizing 
the momentum with respect to kp, the Fermi momentum, the explicit expression for the 
correlation energy per particle in Rydbergs (Ryd.) has the following form, 



E r 



N 167T 5 

3 



d 3 q d 3 p 1 d 3 p 2 



[l-p 2 )e(l- Pl )9(E qi -l)B(E q2 -l) 



q 2 (q + p 1 - p 2 ) 2 (q 2 + q ■ p\ - q ■ p 2 ) 

/oo poo 
dx / dy y 2 {In [1 + H(x, y) + M (x, y) - 21I(x, y)} 
-oo JO 



A-RoPr 2 

-H(x,y) - M(x,y) + 2U{x,y)} 
3 



2na 2 r\ 



/oo poo 
dx / dy y 2 {In [1 + K(x, y) + S(x, y)\ - /C(x, y) - S(x, y)} , (107) 
-oo JO 



where we have set E qi = \q + pi\ and E q2 = \q — p 2 \. 
For convenience, we set 



E c 

— — I t + I 2 + I 3 
N 



(108) 



where I\, I 2 , and Is are the first, the second and the third integrals, respectively. Expressions 
for few of the functions that appear above are given below. The functions JC(x,y), H(x,y) 
and Y,i(x,y) are given by, 

7r(— ix + y l — 1) 
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Figure 10: Ring diagrams generated by II (p) in Eq. (PQ). 




Figure 11: Diagrams that appear in the second In term in Eq.( 
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a. 




b. 

Figure 12: A term of order e 4 that is due to Si. Both representations are equivalent. Initially 
two particles interact with each other with one of them going back to its initial state while 
the other one interacts with a third particle before both returning to their corresponding 
initial states. 



38 



Figure 13: Some of the ladder diagrams due to Ei. 
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Figure 14: Few terms that appear due to S 2 . 
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Figure 15: A term that is both due to Si and E 2 . Here three particles are interacting 
pairwise with two of the three particles exchanging states. 
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Figure 16: Coulomb interactions that must be taken into account in inhomogeneous media. 
These diagrams vanish in the homogeneous case. 
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T-C(x,y) = -e 2 A(x,y)D(x,y) 



(110) 



or more explicitly, we have 
2ar 



H(x,y) 



iry' 



l + -^(x 2 + 4y 2 (l 
8y 6 V 



y 



)(i 



In 



x ( 
arctan 

2y V 



arctan 



2/ 

2 "'" ' 2 

2/(1-1) 



:i) 2 +y 2 (i + i) 2 
:i) 2 +2/ 2 (i-i) 2 



an) 



and 



= 4 B / d 3 k 1 d 3 k 2 g(x,y,k 1 ,k 2 ) 



where 

g(x,y,k h k 2 ) 



(112) 



1 1 (0(1 - fei) - 0(1 - fcia)) (9(h - E') - 0(1 - k 2 )) 

X ■ — 



\y- \y- k 2 \ 2 



-ix 



+ y 2 - l)(ix - k\ - k\ + E' - 1) 



(x) is the step function. E', g(y) and a are given by 

E'= \h + k 2 -y\ 2 , 



(113) 



(114) 



and 



f(?/) 



i + i^m 

2y 



1 + 2/ 



1-2/ 



(115) 



The Gell-Mann-Brueckner term can be isolated from the full expression for the correlation, 
Eq.( |107| ), by rewriting the second integral, I 2 , in the following way: 



jring 



_|_ jexch _|_ jseRing 
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(116) 



where 



jring 



rexch 



3 1 

4-7T r 2 a 2 
3 1 
4ir r 2 a 2 



+ 00 /"OO 

dx J dy y 2 {In [1 + TC(x, y)] — H(x, y)} 
o 

M{x,y) 



-oo 

OO POO 

dx / dyy l \ In 

oo Jo 



1 + H(x,y) 



M(x,y) } . 



(117) 
(118) 



jseRmg - g j g Q f J 2 _ jrmg _ jexch _ rjv^g £ erm j g i nc [eed f u ][ expression for 

the RPA term jjE^ PA . By taking the limit r s — > and y < 1, it reduces to the Gell-Mann- 



Brueckner result 117 



N 



0.0622 In r.-0. 142 + ... 



as we show next. Going back to the expression for T~C(p), Eq.( |99[), we have 



(119) 



jring 



1 /" (i 3 /c f duj 



(2tt)* 



— Ml + H(.k,u))-H(k,v)] 
in 



(120) 



with H(k,u) defined to be 
H(k,u) = 2e 

Now we set, 



2 An f d 3 p (W+kf-\p 



k 2 J (2tt) 3 u2+ (±(p + k)2 _ip2 



2p 2 47T 

w) = — p p + W )) 



(121) 



(122) 



Terms other than J"™ 9 in Eq.( |107| ) provide new correlations to the RPA-term. The first 
term in Eq.( |107| ), as we mentioned above, is the second order exchange term and has been 



evaluated exactly in Ref . ||26|| , 



E^ ch = hn2- JLc(3) (123) 
- 0.04836 Ryd. 
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To compare our expansion to others we calculate the term that is equivalent to the RPA 

calculation. We start by evaluating the integral Hi. After integrating the angular variables, 
we have 



Hi{u, k) = ix / — - — In 
Jo k 



u 2 + {pk + \k 2 ) 2 



oo 2 + {-pk + \k 2 ) 2 
This integral can be easily performed over p. We get after that 



(124) 
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(125) 



The second integral H 2 in Ti(uj, k) is seen to be simply obtained from Hi(cu, k) by replacing 
k by —k. Therefore we have 



H(u,k) = 2Hi(u,k) 



(126) 



If we set y = and x = j-, we get 



H(x,y) = 



2ar s | , , x 2 + 4y 2 (l-^) in 
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arctan 
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(127) 



where r s = . Hence the ring diagrams' contribution to the energy is given by 



E RPA = 3iV 



poo poo 

/ t/ 2 dy / dx [ln[ 1 + H(x, y) ] - H(x, y)\ , (128) 
Jo Jo 



where N is the total number of electrons. Now, we notice that if we make the following 
substitutions: 
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y = y' 

x = 2yx' 



(129) 



we recover exactly the same expression for the ring diagrams as obtained by Bishop and 

Luhrmann in Ref.[25[| using a totally different expansion. This helps to validate our original 
expansion in h and later in iterating the equations of motion in terms of e 2 . We stress that 
this expansion is valid for both small and large r s . 

From the above analysis, we see explicitly that the method of effective action amounts 
to including another infinite set of diagrams besides the usual ring diagrams. One subset of 
the diagrams added is the ring diagrams that allow exchange in them, Fig. [5| On physical 
grounds these are expected to be the next important ones that must be summed up. It is 
also easily seen from the above that all second order diagrams are included in the expansion 
with the right symmetry factors. This shows that our original expansion in ft is indeed 
meaningful. One last thing to note about the diagrams in Fig. |ll] is that they include the 
"anomalous" diagrams |Z7). We have already dealt with these diagrams in the previous 
section. From a nonzero temperature calculation, we were able to show that these diagrams 
give a zero contribution at zero temperature simply because they violate Fermi statistics. 
From Eq.(107), the contribution of these diagrams involves integrating In [1 + K.(x, y)] — 
JC(x, y) over all x in the complex plane. Since K.(x, y) has a simple pole, Eq.(109), the above 
integrand ends up having poles of order two and higher and hence their zero contribution at 
zero temperature. However these diagrams become essential at non-zero temperature where 
the above constraint of Fermi statistic is no longer an issue. In the following section we 
calculate the contribution of the function A4(x, y) to the correlation energy. 
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7 The Inclusion of Second Order Exchange Effects 
in Ring Diagrams 



In this section, we give the contribution of diagrams like those shown in Fig. 17 to the 
correlation energy at zero temperature. Hence we go beyond RPA in this case. We will 
show that our method provides excellent agreement with fully numerical calculations. We 



also compare our results to those found in Ref.p5| where a coupled cluster formalism has 
been used to get the correlation energy. Our final results clearly show that our method is 
more transparent than other previously used methods. To get to these final results we had 
to make approximations along the way. We use two different approximations: the Hubbard 
approximation in Ref.p8| and the Bishop-Luhrmann approximation, Ref. ||25||. We have found 



that the former applies well to high values of r s while the latter applies well to low values of 
r s . From Eq.(118), this amounts to finding J^ xch and the second order exchange term that 
has already been calculated by Onsager et al. in Ref. [pfifl . Now we show how to calculate 
jexch r-Q^ calculation is straightforward but special care must be exercised when it comes 
to numerically evaluating the final result. First we let 



where 



dx I y dy In 
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1 + H{x,y) 



(130) 



M(x,y) 



16a 2 r? 1 



d 3 k\d 3 k 2 



1 
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(131) 



and where the region of integration V is given by: 



r = 



{{k 1 ,k 2 )k 1 < 1, k 2 < 1, \y+h\ > 1, \y- k 2 \ > lj 



(132) 



It is obvious from the above that we are faced with a daunting task of having to deal with a 
9-D integral inside a logarithmic function which itself has to be integrated over normalized 
energy x and normalized momentum y. Most of the complications are related to the angle 
integrals which can not be separated. To be able to make some progress we have to make an 
approximation, i.e., we either assume that \ki + y\ ~ 1 and average over the angle between 
them, i.e., the Hubbard approximation or we use a more sophisticated approximation like 
the one proposed by Bishop and Luhrmann (BL). Hubbard's approximation amounts to the 
following: 



1 



|fci - h + y\ 2 y 2 + 1 



(133) 
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On the other hand, the BL-approximation is more complicated and the reader is referred to 
their appendix in Ref. |25j for a discussion of their approximation. The two approximations 



are almost equivalent for large values of momentum y, i.e. y > 2 in units of Fermi momentum. 
For values of y less than 2, both approximations are quite different ( see Fig. pl| ). The 
Hubbard approximation seems to apply well for moderately large values of r s while the BL- 
approximation applies well for low values of r s . In fact the BL-approximation incorporates 
the Fermi statistics of the interacting particles where it matters most, i.e. in high density 
situations, and this is one reason why we found it to be a good approximation to the exchange 
Coulomb interaction for low r s . The Hubbard approximation should be expected to be a 
good approximation to particles near the Fermi surface and where statistics is not an issue. 
The case of low density should apply well to this case. Given either approximation we end 
up with an expression for A4(x,p) of the form: 

16a 2 r s 2 1 r/ __ N f j3; - x 2 + \(y 2 + 2h ■ y)(y 2 - 2k 2 ■ y) 



( 27r ) y J? (x 2 + \{y 2 + 2h ■ y) 2 j (x 2 + \{y 2 - 2k 2 ■ y^j 

where the function f(y) is our approximation to the exchange term in M. (x, y) and is assumed 
known. 

The integrations over ki and k 2 are easily found if we choose to write both vectors in 
cylindrical coordinates with y along the z-axis, i.e., 

h = (pi,Zi,8i), i= 1,2. (135) 
Hence the integration over k\ becomes: 



2- 



d k\ = / d6i I pi dpidzi , (136) 



where 



T' = {y<2, -l <Zl< l- y , (l-fa + y) 2 )* < Pl <(l-z 1 2 ) h } 
\j{y<2, l-y< Zl <l, 0< Pl <(l-z 1 2 )"} 

\j{y>2, -1<Z!<1, 0< P i < {l-z^Y (137) 

with an equivalent expression for the region of integration over k 2 . The integrations are now 
easily carried out. The expression that we get for A4(x, y) is naturally expressed in terms of 
Y instead of x. So in the following x refers to We only quote the final expression here: 
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|iz£ ( arcta n(^) _ arctan^)) jj) 0(2/ - 2) . (139) 



Given this expression for Ai(x, y), we can easily write the full expression for the correlation 
due to exchange effects in ring diagrams. It is important to note that we allow for more 
than one exchange to take place in each ring diagram. Hence our original integral I, can be 
now evaluated numerically. So within the above approximation, we were able to reduce our 
calculations to a two-dimensional integral: 



jexch 
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(140) 



The expression for the function f(p) depends on which approximation we choose to use. 
In this formalism it is hard to decide on which one since both give very close answers to 
the second order exchange diagram. We have found that the Hubbard approximation in 
this instance gives 0.041 Ryd. while the BL approximation gives 0.049 Ryd. Clearly the 
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latter is closer to the true value of 0.0484 Ryd. } but we believe that this is not enough to 
decide against using the Hubbard approximation. In fact at low r s we expect that the high 
energy part of the exchange energy between particles to dominate and this explains why 
we found the BL-approximations to apply well in this range. At high r s this is not true 
anymore and in this instance we expect the BL- approximation to overestimate this exchange 
of energy between the particles. This is in fact what we have found. This casts doubts that 
the BL-approximation is always better than the Hubbard approximation. As we said above 
the second order exchange term is not a good test since it is density independent. This 
explains why the Hubbard approximation gives better results for higher values of r s . Below 
we give some numerical answers to compare with those obtained using quantum Monte Carlo 
methods (QMC) fl29| . Our results are obtained using both approximations. For the second 
order exchange diagram we have used the true value of 0.0484 Ryd. 

Clearly our results, Fig.[l9], agree well with previous results. For large values of r s , i.e., r s 
~ 100, our results do not compare well to the QMC results. However this is expected since 
we need to include corrections due to the ladder diagrams which are believed to be important 
at low densities. This concludes our calculations. We have shown the effectiveness of an 
effective action approach to electronic systems. In the next section we shall state the viability 
of using this formalism in the nonhonogeneous case and show how it is related to density 
functional theory. 
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Figure 19: Comparison of final results for the correlation energy (in Ryd.) of a homogeneous 
Fermi gas. 
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8 Conclusion: The Nonhomogeneous Electron System 



In this last section, we would like to go back to the inhomogeneous case and show how our 
method gives some of the results known in DFT || ||. We start by briefly reviewing the 
basic ideas behind DFT. The Hamiltonian H is written as the sum of three terms: T, V 
and U. T is the kinetic energy term, V is the external potential term and U is the Coulomb 
energy term. So we have 

H = T + V + U , (141) 
T = \ J W> f (r)W>(r)d 3 r , (142) 

V = V(r)^(r)i(j(r)d 3 r , (143) 



u= ir f(rWp^^ dr dr , 

2 J \r — r'\ 

If \1/ is the ground state, then the density n(r) is given by 

n(r) = (y\^(r)4>(r)\ty . (145) 

The first important fact to note is that V(r) is a unique functional of n(r), assuming that 
there is a unique ground state for the system. Now let the functionals F[n(r)] and Ey[n(r)] 
be defined as follows: 

F[n{r)} = {V\T+U\V) , (146) 

E v [n{r)} = J V(r)n(r)d 3 r+ F[n(r)) . (147) 

Then, if we assume that there is a one-to-one correspondence between n(r) and V(r), it can 
be shown that 

5Ev\n] . 

wi»- = ' (148) 
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with 

n(r) d 3 r = N 



This result is called the Hohenberg-Kohn theorem. Next, the Coulomb energy is isolated 
from the functional F[n] by introducing a new functional G [n\. 

Fin] = - [ drd7 U ^ U ^ + G\n] , (149) 
z J \r — r'| 

G[n] = \ [ dr V r V r , m(r, r')\ r=r , + \ [ ^^-dr dr', (150) 
Z J Z J \r r 



and 



C 2 (r, r') = n 2 (r, r'\ r, r') — n(r)n(r') . (151) 



Here n\{r, r') is the one-particle density matrix and C^r, r') is a correlation function defined 
in terms of the one- and two-particle density matrices. DFT calculations are essentially 
centered around finding good approximations to the functional G[n]. This is usually done 
by postulating that there is a virtual system of free electrons in an external potential with 
exactly the same density as the interacting system. The energy is found by finding the 
eigenfunctions that correspond to this external potential self consistently. 
The method that was presented here is in fact similar in many ways to the above ideas 
expressed in DFT. However there are major differences. Our method clearly incorporates 
the Hohenberg-Kohn theorem. From Eq.(35), we have the following 



From Eq.(49), the term 



|SiM = -QM. (152) 
op(x, x) 



d 3 x V (x)p(x, x) , (153) 



can be separated from the energy E gs . Hence if we set Q(x, x) = , we get right away 
the Hohenberg-Kohn result. In fact we now have an explicit expression for Ey[n] within 
perturbation theory. This result also shows the one-to-one correspondence between density 
and external potential as long as there is a unique solution around (J,B,Q) = (0,0,0) for 
the defining equations, Eqs. (33-35). We have solved these equations only approximately. We 
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have made two important approximations. The first corresponds to the number of diagrams 
included in T from the outset. However from the calculations of the homogeneous case, we see 
that it will probably be the case that including only two diagrams in the inhomogeneous case 
will give good results. Including higher order corrections by taking into account diagrams 
like the ones in Fig. |20] is nevertheless straightforward if more accurate results are desired. 
The second approximation we have made is that in solving for p(x,y). We had to linearize 
the equations of motions to be able to solve for p(x, y) iteratively. Also from Eq.(50), we get 
the following equivalent result: 

VV c (r) = 4vm(r) . (154) 



Clearly our method is a generalization of the Thomas-Fermi method where n(r) in the last 
equation is the true density of the system. A very important difference from DFT is that 
within this method we have a systematic scheme for calculating the functional G[n], Eq.(49). 
From the above analysis in the homogeneous case, it is obvious that the most important 
contribution to correlation at zero temperature comes from calculating the following term: 

TrUX = U to*'* (155) 

where 

X = A- X C . (156) 

To calculate the determinant of the above operators, we need a basis of wavefunctions. In 
practical computations, it is more advantageous to use a finite basis set. If we choose as our 
basis set the wavefunctions {4>i}f =1 with eigenvalues ^ such that 

_I V 2 + V{x)j <j>i(x) = e^ix), (157) 

then the wavefunctions, ipj, of the full interacting system of N electrons can be written as a 
linear combination of the above wavefunctions: 

A' 



ipj(x,t) = ^djifaWe e *'. (158) 



i=l 



The coefficients aji are found self consistently by solving Eq.(56). To solve for X(x,y), we 
solve Eq.(53) self-consistently. Clearly our method is computationally more intensive than 
DFT. However here we are on more secure grounds and our approximations are very well 
controlled to any desired accuracy. Plans to apply this method to simple atoms are underway. 
We hope to report on these developments elsewhere. Finally, we would like to point out that 
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Figure 20: Diagrams that appear when we include terms of order h 6 in T. 
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this method applies equally well to excited states other than the ground state. The Hartree 
field of the excited state is found again by using Eq.(33). If tpo(x) is the ground state solution 
then tpi(x) = Lpc(x) + A(x) is the excited state solution iff 



6T 



(159) 



8ip{x) 



<P=<pi 



0. 



This implies that 




5 2 T 



A(y) = 0. 



(160) 



6<p(x)6<p(y) 



In the above A(x) is assumed small compared to tpo{x). This latter equation enables us to 
solve for A(x) and hence the corresponding p(x, y). Clearly effective actions methods are in 
general superior to conventional methods as is well known in field theory. We expect this 
method to apply equally well to atomic and molecular systems as we have shown to be the 
case for homogeneous systems. 
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